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ON THE MODELING OF SHELLS IN MULTIBODY DYNAMICS 

OLIVIER A. BAUCHAU*, JOU-YOUNG CHOI+, AND CARLO L. BOTTASSO* 


Abstract. Energy preserving/decaying schemes are presented for the simulation of the nonlinear multi¬ 
body systems involving shell components. The proposed schemes are designed to meet four specific require¬ 
ments: unconditional nonlinear stability of the scheme, a rigorous treatment of both geometric and material 
nonlinearities, exact satisfaction of the constraints, and the presence of high frequency numerical dissipa¬ 
tion. The kinematic nonlinearities associated with arbitrarily large displacements and rotations of shells are 
treated in a rigorous manner, and the material nonlinearities can be handled when the constitutive laws stem 
from the existence of a strain energy density function. The efficiency and robustness of the proposed ap¬ 
proach is illustrated with specific numerical examples that also demonstrate the need for integration schemes 
possessing high frequency numerical dissipation. 

Key words, multibody dynamics, geometrically exact shell, time integration, energy decaying scheme 

Subject classification. Applied and Numerical Mathematics 

1. Introduction and Motivation. This work is concerned with the numerical simulation of geomet¬ 
rically exact shell models within the context of multibody system dynamics. While the partial differential 
equations that govern shell problems are well known and presented in numerous textbooks, their numerical 
treatment is still the subject of active research. Indeed, numerical analysis tools for partial differential equa¬ 
tions have significantly changed in recent years. In the past, general purpose discretization methods were 
developed, with emphasis on robustness, performance, and accuracy. These methods aimed at solving vast 
classes of problems such as ordinary differential equations, differential/algebraic equations, or hyperbolic 
conservation laws. 

This approach is now changing. The differential equations that govern many problems in mathematical 
physics possess qualitative and structural characteristics that can be determined by studying their geometry. 
Classical examples of such characteristics are the invariants associated with Hamiltonian systems, the sym- 
plectic structure of the governing equations, or symmetries and attractors. There is increasing evidence that 
numerical methods that correctly recover the qualitative features of the underlying differential equations are 
often endowed with superior computational performance, greater robustness and improved accuracy. 

This new paradigm has resulted in the development of a new mathematical discipline, called geometric 
integration , whose aim can be summarized as: “the purpose of computing is insight, not numbers”, to quote 
Hamming [17]. Geometric integration theory is now becoming a bridge that links the work of pure, applied 
and computational mathematicians, and it is usually expressed using the terminology and formalism of 
differential geometry and Lie group theory. 

In reality, engineers have used geometric integration for a number years. Simo and his co-workers were 
among the first to develop special integration procedures for nonlinear structural dynamics. They analyzed 
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the problem of the dynamics of rigid bodies [24], nonlinear elasto-dynamics [21], geometrically exact shells [22] 
and geometrically exact beams [23]. In all cases, the idea was to design algorithms that ensure the discrete 
preservation of the total mechanical energy of the system, therefore obtaining unconditionally stable schemes 
in the nonlinear regime. These energy methods, sometimes equipped with additional conservation properties, 
such as the conservation of momenta, are becoming increasingly popular in multibody dynamics, paralleling, 
although unknowingly, the mathematical developments of geometric integration theory. 

However, increasing evidence points toward the fact that geometric integration is not sufficient, per se, 
to obtain robust integration schemes. While these schemes perform well for problems with a small number of 
degrees of freedom or those featuring a “smooth” dynamic response, they tend to be quite unsatisfactory when 
applied to the complex simulations encountered in many engineering applications [6]. In fact, the predicted 
time histories of internal forces and velocities can present a significant high frequency content. Furthermore, 
the presence of these high frequency oscillations hinders the convergence process for the solution of the 
nonlinear equations of motion. The selection of a smaller time step does not necessarily help, as a smaller 
time step allows even higher frequency oscillations to be present in the response. These oscillations are 
particularly violent in multibody dynamics simulations because these systems are rather stiff due to the 
presence of numerous algebraic constraints, while the nonlinearities of the system provide a mechanism to 
transfer energy from the low to the high frequency modes. These difficulties can easily go unnoticed in 
many test cases, but will create major problems for large scale simulations. Consequently, the presence of 
high frequency numerical dissipation appears to be an indispensable feature of robust time integrators for 
multibody systems. Such feature can be added to schemes developed within the framework of geometric 
integrators, as it will be shown in this paper. High frequency numerical dissipation does not alter the physics 
of the problem, since the high frequency content of the response filtered out by the algorithm is itself an 
artifact of the spatial discretization process and contains no information about the physical behavior of 
the system. The need for high frequency dissipation in large finite element models has been recognized for 
many years in various disciplines, from structural dynamics to fluid mechanics. The most widely used time 
integrators for finite element analysis are high frequency dissipative [19, 18]. 

This paper focuses on the development of a geometric integrator for shell structures that preserves 
important qualitative features of the underlying equations, and is equipped with high frequency numerical 
dissipation. The goal of the work is to obtain schemes presenting improved robustness and reliability over 
standard “black-box” integrators. In order to achieve this goal, the specific features of the equations governing 
nonlinear flexible multibody systems with shells are reviewed. 

First, the governing equations are characterized by linear and rotational tensorial fields describing kine¬ 
matic (displacements, velocities) and co-kinematic (forces, momenta) quantities. In shells, the rotational 
field describes the evolution of unit director, and is therefore a special family of two-parameter rotations. 
Second, the equations are nonlinear because of large displacements and finite rotations (geometric nonlin¬ 
earities), and possibly because of nonlinear constitutive laws (material nonlinearities). Third, the presence 
of joints imposes different types of kinematic constraints between the various bodies of the system. In this 
work, the Lagrange multipliers technique is used to enforce the constraints, giving the governing equations a 
differential/algebraic nature. Fourth, the equations of motion imply the preservation of a number of dynamic 
invariants, in particular the total mechanical energy, and the total linear and angular momenta. 

Since discrete preservation of energy leads to unconditional nonlinear stability, preservation of this in¬ 
variant is a central focus of this paper. The proposed geometric integration procedure is therefore designed to 
satisfy specific requirements. First, a discretization process is developed that preserves the total mechanical 
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energy of the system at the discrete solution level, as well as the total linear and angular momenta. This 
process is independent of the spatial discretization procedure that is left arbitrary. In the present implemen¬ 
tation, the finite element method is used, and the mixed interpolation of tensorial components [1, 2, 16] is 
implemented to avoid the shear locking problem. Next, the reaction forces associated with the holonomic 
and non-holonomic constraints imposed on the system are discretized in a manner that guarantees the sat¬ 
isfaction of the nonlinear constraint manifold, i.e . the constraint condition will not drift. At the same time, 
the discretization implies the vanishing of the work performed by the forces of constraint at the discrete 
solution level. Consequently, the discrete energy conservation laws proved for the flexible members of the 
system are not upset by the introduction of the constraints. The resulting Energy Preserving (EP) scheme 
is a geometric integrator for multibody systems with shells that provides nonlinear unconditional stability. 
However, it clearly lacks the indispensable high frequency numerical dissipation required to tackle realistic 
engineering problems. 

Using a simple procedure [3, 4] based on the EP scheme, it is possible to derive a new discretization 
that implies a discrete energy decay statement. In the resulting Energy Decaying (ED) scheme, the system 
no longer evolves on the constant energy level set, but is allowed to drift away from it in a controlled 
manner. This concept seems to be new in geometric integration theory, and provides a procedure to obtain 
nonlinear unconditional stability (from the bound on the energy), together with a mechanism for removing 
the undesired high frequencies. The discretization process for the forces of constraint is left unchanged: 
the work they perform vanishes exactly, while the system evolves on the constraint manifold without drifts. 
Therefore ED schemes satisfy all the requirements set forth earlier. 

Related ED schemes for various problems in multibody dynamics were proposed in the literature. Finite 
difference schemes were presented in [8, 9, 5, 13]. Time discontinuous Galerkin approximations of the 
equations of motion written in the symmetric hyperbolic form are used in [7, 3]. In [14], ED schemes are 
cast in the form of Runge-Kutta schemes and related to the basic concepts of geometric integration theory. 
The analysis in ref. [4] shows that slightly different EP and ED discretizations can be developed, usually 
through different treatments or parameterizations of finite rotations. Some discretizations might also present 
additional conservation properties. For instance, some EP and ED schemes also imply the conservation of 
momenta, or are geometrically invariant [10, 15, 4]. These additional features are easily obtained by recasting 
the field equations in fixed pole form [11], a procedure that brings back again the link with Lie groups through 
the use of exponentials and of Cayley’s transform. 

The paper is laid out as follows. The equations of motion for shells are presented in section 2. The energy 
preserving scheme presented in section 3 is then generalized to an energy decaying scheme in section 4. A 
special element, the shell revolute joint, used to connect shells to other components of a multibody system, 
is developed in section 5. Finally, numerical examples are presented in section 6. 

2. Formulation of the Equations of Motion. 

2.1. Shell Kinematics. Consider a shell of thickness h and reference surface area ft, as depicted in 
fig. 2.1. An inertial frame of reference S consisting of three mutually orthogonal unit vectors ii, i 2 , h is 
used. Let be the position vector of an arbitrary point on the reference surface of the shell, and let £ be 
the material coordinate along n, the normal to the reference surface. The position vector r of an arbitrary 
point on the shell in its reference configuration is then 

de ,e,o=L 0 (t l ,?)+c ,e), (2- 1 ) 
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where ^and £ 2 are the material coordinates used to represent the shell reference surface. The coordinates 
£ 1 , £ 2 and ( form a set of curvilinear coordinates that are a natural choice to represent the shell geometry. 
The coordinates £ x and £ 2 are assumed to be lines of curvatures of the shell reference surface. The base 
vectors are then 



< 9 r' 


£i> £2’ £3 


— 


(1 ” £li (1 ) —2 > — 


( 2 . 2 ) 


where R\ and R 2 are the principal radii of curvature, a a — and the notation (•),<* is used to denote a 
derivative with respect to £ a . It is convenient to introduce a set of three mutually orthogonal unit vectors 
at the shell reference surface (i.e. at £ = 0) 


ei = 


vsir 


e 0 = 


^2 . 
\/tt22 ’ 


e 3 “ n, 


(2.3) 


where a aa = a a ■ a a . 

According to the classical inextensible director model, it is assumed that the material line initially normal 
to the reference surface of the shell remains a straight line and suffers no extension. Therefore, the position 
vector of a material point of the shell can be written as 


0 




(2.4) 


where u(£ x ,£ 2 ) is the reference surface displacement vector. The base vectors at the shell reference surface 
in the deformed configuration are 


G = [G l5 G 2 , G s \ — 





(2.5) 


Introducing the position vector, eq. (2.4), then yields 


G = 


[a g 2 

— . ~~, ~~ , O’ 3 

.V a H V a22 


= E + (H, 


( 2 . 6 ) 
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where 


E — [E x , E_ 2 ) £ 3 ] — 


£1 + 


y/an' 


£2 + 


-=Bz, £3 

V^2 


= 


£3,1 £3,2 q 

\^n’ v^2 5 J 


(2.7) 


Note that £ 3 (^ 1 ,^ 2 ) is a unit vector, whereas and E 2 are not unit vectors, nor are they orthogonal to 
£ 3 , as axial and transverse shearing strains develop during deformation. 


2.2. Equations of Motion. The Green-Lagrange strain tensor e is defined as 

e = ^(G T G - g T g). (2.8) 

The strain tensor e is defined in the curvilinear coordinate system defined by coordinates ^ 1 ,^ 2 and £. 
However, it is more convenient to work with strain tensor e defined in the locally rectangular system defined 
by triad e x , e 2 , eg, see eqs. ( 2 . 3 ). For shallow shells ( i.e . C/Ri 1 and (/R 2 1) undergoing large 

displacements and rotations but small strains (all strain components are assumed to be small compared to 
unity), the strain-displacement relationships can be written as 


where 


e = i [E t E - I + C {E t H + H t E + k)] , 


K> = 


l/Rr 0 0 

0 1/R 2 0 

0 0 0 


(2.9) 


( 2 . 10 ) 


It is clear that the strains can be expressed in terms of five parameters: the three components of the 
displacement field u (through E x and E 2 ) and the two parameters defining the orientation of the unit 
director E s . Virtual changes in the strain energy of the structure are given by 

SV= [ j SV dCdft = [ [ Se r dCdH, ( 2 . 11 ) 

JQ Jh Jn Jh 

where SV is the virtual strain energy density, and r the second Piola-Kirchhoff stress tensor. Introducing 
the strains, eq. ( 2 . 9 ), and taking into account the symmetry of the stress tensor then yields 


SV = SE ■ {E + C H)r + SH • (Er. 


( 2 . 12 ) 


The existence of a strain energy density function V is postulated here, hence the constitutive laws are of the 
form r = dVjde. 

The velocity vector of material point P of the shell is obtained by differentiating the position vector, 
eq. ( 2 . 4 ), with respect to time, to find v = u + (£ 3 * The kinetic energy of the system is now 

K=[ [ K dCdfi = \ ( [ pv-v dCdfi, (2.13) 

Jn Jh 2 Jn Jn 

where K is the kinetic energy density. Introducing the velocity vector then yields 

K=^p{u + CE 3 )-{u + CE^). (2.14) 
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Hamilton’s Principle now writes 


[ ! [ [ (SK — SV ) d<dfld t 
JJ Q Jh 

= J*' J J [p(6u + (6E ?> ) ■ (u + CI 3 ) +SE-(E + C H)t + SH ■ QEt] d(dndt 


= 0. (2.15) 


Integrating through the thickness of the shell, we get 


J*' J {<5« • [k - (N hl + jV 2 , 2 )] + 5E - [g f - (M 1;1 + M 2 . 2 ) + Es ]} 


dCldt = 0 . 


In this expression, h = mu + and g = + /*£ 3 are the linear and angular momentum vectors of 

the shell, respectively; the mass coefficients are defined as m = f h p d£, s* = f h p( d(, I* = f h p( 2 d(. The 
spatial in-plane forces are N a = (EjV* + HM.* a )/y/a W, the out-of-plane forces N 3 = and the bending 
moments M a = (EM* a )/^a^. The convected forces are iV* = [iV^,iV*,iV*] = f h r dC, and the bending 
moments M* = [M*, Mo, M 3 ] = f h tQ d£. 

3. Energy Preserving Scheme. Discrete equations of motion that imply conservation laws for the 
total mechanical energy, linear momentum and angular momentum of the system will now be developed. 
Times U and tf denote the initial and final times for a time step, respectively, and the subscripts (•)* and (•)/ 
indicate quantities at U and tf, respectively. Furthermore, the subscript («) m is used to denote mid-point 
average quantities defined as 

(•)m = | [(•)/ + (•)<]- (3- 1 ) 

The following matrix identity will be used extensively 

AjB f - AjBi = (A f - Ai) T B m + A T m (B f - B { ). (3.2) 

Hamilton’s principle, eq. (2.15), is now approximated in time in the the following manner 

Sj\p [<*/ - «<> + «&/ - a.)] • + < 

+(j Ef — Ei) • ( E m + C H m )T a + (Hf — Hi) ■ d£dfi = 0. (3.3) 

The change in strain components from ti to tf is evaluated with the help of identity (3.2) to find 

e f - a = \ [(E f - Ei) T (E m + (H m ) + (E m + C H m f(E f - E { ) 

+(H f - Hi) T (E m + C El(H f - Hi)] . (3-4) 

Over one time step, the strain components can be approximated as e(r}) — e m 4- 77 (e/ — e«)/2, where 77 = 
2(t - t m )/ At is the non-dimensional time. If the strain energy density function V is viewed as a function of 
the scalar variable 77, the mean value theorem then implies the existence of a fj G [— 1 , 1 ] such that 

Vf = Vi + - 7 T— —2 = V{ + r a ♦ (e/ — e^. (3.5) 

J de n dg 
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This relationship defines the average second Piola-Kirchhoff stress tensor, r a = dV/de\ Combining this 
result with eq. (3.4) then leads to 

{E f - Ei) ■ (E m + C H m )r a + (H f - Hi) • C = (e/ - «) -T a = V f - V u (3.6) 


where the symmetry of the stress tensor was taken into account. For linear constitutive laws of the form 
r — C* e, where C* is the stiffness matrix, the average stress tensor simply becomes r a — C* e m . 

The following configuration updates are now defined 


Uf-U* _ . 

At “ m5 


£ 3 / 


Eo 


At 


= E- 




(3.7) 


Introducing eqs. (3.6) and (3.7) into the approximate expression for Hamilton’s principle, eq. (3.3), then 
leads to 



(u m + c MsJ • [(u f - Ui) + c {Ms/ - fij] + (V> - Vi)} dCdQ 

J Q J h {%(u f + CEsf) • (u f + C Is/) - f (Ui + CEsi) • (Mi + C £Ui) + (Vf - Vi)} dCdn 


= [ [ [{K f - Ki) + {Vf - Vi)] dCdf2 = 0. 
Jn Jh 


(3.8) 


This result clearly implies the conservation of the total mechanical energy of the system within a step. 
In summary, the approximate form of Hamilton’s principle given by eq. (3.3) leads to a discrete energy 
conservation statement, eq. (3.8), when the configuration updates are chosen according to eqs. (3.7), and the 
average stress according to eq. (3.5). 

Integrating through the thickness of the shell leads to 




Ef ~ hi 


At 


- +iV 2m ,2) 


+ (£ 3 / - Ssi) * 


-+ M.2m$) + —3m | 


At 


dO = 0. (3.9) 


In this expression, the in-plane forces are N_ am = (E m N* aa + H m M-* aa )/ the out-of-plane forces N^ m = 
E m N^ a , and the bending moments M am — (EmM* aa )/y/a aa . The discrete governing equations of motion 
for shells are then 

Ql, - Ql iM„„, + - a !ln ) = (3.u) 

where p are the externally applied loads, and q* the externally applied moments measured in the local 
system. The finite change in director orientation E^f — E^ { was expressed in terms of the two parameter 
incremental rotation vector, see B. 

Invariance of the system Hamiltonian under spatial translations and rotations implies the conservation 
of the linear and angular momenta. The preservation of these invariants is less crucial from a numerical point 
of view, since it does not lead to non-linear notions of stability as in the case of the energy. However, is can 
be of some interest to note that momenta are indeed preserved at the discrete solution level by eqs. (3.10), 
(3.11) and (3.7). The complete proof of this assertion is reported in ref. [12], and it is based on projecting 
the discrete equations of dynamic equilibrium onto suitable test functions and integrating over the shell. 
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Finally, it is important to point out that the particular spatial discretization adopted plays no role in 
the proofs leading to the properties of discrete energy and momentum conservation. Therefore, any spatial 
discretization of the discrete equations of motion will inherit these properties, when the configuration updates 
are chosen according to eqs. (3.7), and the average stress according to eq. (3.5). 

4. Energy Decaying Scheme. In this section we derive a scheme associated to a discrete law of 
energy decay, using the EP scheme as a basic building block. The resulting integrator for shells features high 
frequency numerical damping, overcoming the difficulties and lack of robustness of EP methods. 

First, an additional state is introduced at time tj = \im € -^o(U + e), and the subscript (*)j is used to 
denote quantities at this time. The following averages are now defined 

(•), = ^ [(•)/ + (*)j]; (*)h = 2 [(*)i + (•)<] • 

The ED scheme proceeds from the initial to the final time by means of two coupled steps: one step from tj 
to tf, the other from t, to tj. The time-discrete equations of dynamic equilibrium are 


Tt - (^l9,l + ^2 9 ,2) - Pg ; 

Qm - <% (Ml,. 1 + M 29 ,2 - Z 3g ) = £*; 


(4,2) 


° + 3 [(=^19,1 + ^29,2) “ (i^lp.l + —2p, 2 )] - P h \ 

qi g± xf L + \[q t 9 (Mi,,i+m 2 9, 2 - K 3g ) - qi (Mip.1 + m 2p<2 - m P )] = «;• 


(4.3) 


At 3 

Q j-Si 

At ' 3 

The configuration update relationships are given as 

Uf = Wj + At (u f + Uj)/ 2, Uj = Ui - At [u f -Uj- a(uj - «,)] /6; 

E$f = £Ui + At (£ 3 / + 4j)/2) £ 3 j = On ~ At [4/ — 4» — a(E 3 j — 4<)] /6, 

where a is a tuning parameter that controls the amount of numerical dissipation provided by the scheme, 
while the forces N ap and moments M ap are given by 

N ap = N ah + a (N aj - N_ai)/2] M ap = JkU + a (M aj - M ai )/2. (4.5) 


(4.4) 


Using developments similar to those exposed for the EP scheme, it can be easily shown that the proposed 
discrete equations imply 


(. K f + V» - (Ki + Vi) + a c 2 = 0. 


(4.6) 


c 2 is a positive quantity given by 

c 2 = /J [m II u II • II u II +2s* II u II • II 4 II +r II 4II • II4 ll] M 

+ j\ II e II C7* IS e II do > 0, (4.7) 

where || • ||= (*)j — (•)* is the jump between and tj. This result implies the decay of the total mechanical 
energy over one step of the algorithm, (Kj + Vf) < (K{ + V*). The parameter a can be used for controlling 
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Fig. 5.1. Configuration of the shell revolute joint in the reference and deformed configurations. 

the amount of energy that is dissipated within the step. It should be noted that the property of preservation 
of momentum observed in the EP case is lost in the ED algorithm. 

It is interesting to characterize the ED scheme with the classical tools of linear analysis. In fact, if the 
scheme is applied to a single degree of freedom linear oscillator model problem, the asymptotic value of the 
spectral radius of the amplification matrix, />oo, is found to be poo = (1 — ot )/(1 + a). For a = 1, poo = and 
asymptotic annihilation is achieved. If a = 0, poo = 1? and i n view of eq. (4.6), energy is exactly preserved. 
Hence, the ED scheme is in fact a family of schemes with a single tuning parameter, a, that controls the 
amount of high frequency numerical dissipation; both asymptotic annihilation or exact energy preservation 
can be achieved with the same scheme by using a = 1 or 0, respectively. 

5. The Shell Revolute Joint. When modeling multibody systems with shell components, connections 
between shell elements and other elements of the model must be carefully treated. Indeed, beam or joint 
nodes involve six degrees of freedom, three displacements and three rotations. On the other hand, shell nodes 
involve five degrees of freedom, three displacements and two rotations that determine the orientation of the 
director E 3 . A typical situation is shown in fig. 5.1 that depicts the connection between a beam and a shell. 
If the beam were clamped to the shell, the beam would apply bending moments and a “drilling moment” 
at the connection point. The shell cannot sustain such drilling moment since it presents no stiffness about 
an axis normal to its reference surface. Hence, the connection between the beam and the shell must be 
done through a shell revolute joint , i.e. a revolute joint with its axis of rotation perpendicular to the shell 
reference surface. 

Consider a beam and a shell denoted with superscripts (.)* and (.)*, respectively, linked together by a 
shell revolute joint, as depicted in fig. 5.1. In the reference configuration, the normal to the shell is n and the 
triad e 10 , e 20 , and = n is attached to the beam at the connection point. In the deformed configuration, 
no relative displacements are allowed and the beam attached triad rotates to e 1} e 2 , and e 3 = E 3 . This 
condition implies the orthogonality of E^ to both e x and e 2 . The kinematic constraints associated with a 
shell revolute joint are 

C a = E 3 -e a =Q, (5.1) 

where a = 1,2. Of course, at the connection point, the displacements of the beam and shell are identical; this 
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L= lm 



Fig. 6.1. Configuration of the lateral buckling problem. 


constraint is readily enforced within the framework of finite element formulations by Boolean identification 
of the corresponding degrees of freedom. Holonomic constraints are enforced by the addition of a constraint 
potential A C, where A is the Lagrange multiplier. The details of the procedure used to enforce constraints 
in such a manner that the work done by the constraint forces vanishes exactly can be found in refs. [3, 5]. 

6. Numerical Examples. All the examples described in this section will be treated with the proposed 
ED family of schemes corresponding to values of the tuning parameter a 6 [0,1], Although any value of a 
within this range can be used, the examples described here will contrast the two extreme choices. For a = 1 
(p^ = 0), asymptotic annihilation is obtained, and this will be called the ED scheme. On the other hand, 
for a — 0 (poo — exact energy preservation is achieved, and this will be called the EP scheme. 

6.1. Lateral Buckling of a Thin Plate. Consider a thin cantilevered plate acted upon by a crank 
and link mechanism, as depicted in fig. 6.1. The plate is of length L = lm, width h = 0.08 m and thickness 
t = 2 mm. It is clamped along edge AB and reinforced along edge CD by a beam with a square cross-section 
of side a = 4 mm. At point C , the beam connects to a crank and link mechanism. The crank of length 
L c = 0.1 m is attached to the ground at point G , and the link is of length L c = 0.5 m. The ground, 
crank, and link are connected together by means of revolute joints, whereas the beam and link are connected 
through a spherical joint. All components are made of steel with the following properties: Young’s modulus 
E — 210 GPa, Poisson ratio v = 0.25 and density p — 7870 kg/m 3 . The crank rotates at constant angular 
velocity Q = 1 rad/s, and the system is simulated for a period of 2i r s, corresponding to a complete revolution 
of the crank. 

The system is modeled first using a geometrically exact beam element, then using the shell element 
described in this paper. The crank and link are modeled by rigid bodies. For the beam model, three four- 
noded, geometrically exact beam elements were used, whereas for the shell model a 6 x 2 grid of quadratic 
elements was used. The simulation used the ED version of the proposed algorithm (po© = 0.0) with a 
constant time step At = 10 _O3 s. 

As the crank rotates, the plate deflects downwards then snaps laterally when its buckling load is reached. 
In the post-buckling regime, the plate becomes significantly softer in bending due to its large twisting allowed 
by the spherical joint. These features are illustrated in fig. 6.2 that depicts the post-buckling of the plate at 
various times. The plate mid-span lateral deflections obtained for both beam and shell models are shown in 
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Fig. 6.2. System configurations at various time instants during the simulation. 



Fig. 6.3. Time history of the mid-span lateral displacement of the plate. Shell model (ED scheme): solid line; shell model 
(EP scheme): dashed line; beam model: dash-dotted line (results shifted downwards 0.1 m for clarity). 
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Fig. 6.4. Time history of the mid-span lateral displacement of the plate. Shell model (ED scheme): solid line; shell model 
(EP scheme): dashed line (the sign of the displacement was changed to ease the comparison); beam model: dash-dotted line. 
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CRANK ROTATION [rev] 


Fig. 6.5. Time history of the shell mid-span in-plane shear force F 12 ■ ED scheme: solid line; EP scheme: dashed line 
(results shifted downwards 5 10 4 N/ra for clarity). 


fig. 6.3, and good agreement between the models is observed. In both cases, elastic vibrations are superposed 
onto the overall motion imparted by the crank. An enlarged view of the onset of buckling for both models is 
shown in fig. 6.4; the sudden appearance of lateral deflections and subsequent oscillations are observed. This 
figure also shows that the agreement between the beam and shell models is qualitative, not quantitative, as 
should be expected. 

Prior to buckling, the plate resists the bending loads applied by the driving mechanism with very little 
deformations, and high shear forces build up in the plate. When buckling occurs, twisting of the plate 
renders it much softer in the vertical direction, offering little resistance to crank motion. Figs. 6.5 and 6.6 
show the mid-span in-plane shear force F\% and crank driving torque, respectively. Note the linear increase 
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Fig. 6.6. Time history of the crank driving torque. ED scheme: top figure; EP scheme: bottom figure. 

of the driving torque up to Q « 60 Nm, followed by an abrupt drop at buckling. This jump excites the 
vibratory modes of the system. 

Next, the same problem was simulated using the EP schemes, Le. (poo = 1-0). Lateral displacements, 
see figs. 6.3 and 6.4, are found to be in good agreement with the ED predictions. The change in sign for 
the lateral deflection is immaterial since the plate can buckle in either direction. Figs. 6.5 and 6.6 compare 
the predictions of the EP and ED schemes for mid-span in-plane shear force F 12 and crank driving torque, 
respectively. The EP predictions show high frequency oscillations with amplitudes an order of magnitude 
larger than those predicted by the ED scheme. This numerical noise completely obscures the results of the 
computation. 

6.2. Snap-Through of a Cylindrical Shell. A crank and link mechanism is used to drive a cylindrical 
shell through an unstable, snap-through configuration. The system geometry is depicted in fig. 6.7. The shell 
consists of a 60 degree sector of a cylinder of height h = 2.5 m, radius R = 5 m and thickness t = 8 mm. The 
two straight edges of the shell are simply supported, whereas the other two are free. The shell is reinforced 
along line BE by a beam with a square cross-section of side a = 20 mm. At point E, the beam connects 
to a crank and link mechanism. The crank of length L c = 1.5 m is attached to the ground at point G 
located 5 m below point E. The beam, link, crank, and ground are connected together by means of revolute 
joints. The crank is modeled as a rigid body, while the link is a beam with a square cross-section of side 
s — 40 mm. The shell and reinforcing beam are made of aluminum: Young’s modulus E — 73 GPa, Poisson 
ratio v — 0.30 and density p = 2700 kg/m 3 . The lever is made of steel: Young’s modulus E — 210 GPa, 
Poisson ratio v — 0.30 and density p = 7800 kg/m 3 . 

The crank rotates with the following schedule 

(l-cos27rt/T)/2 t <T/2 (61 x 

t > T/2 ’ 1 ’ 

where T — 3 s. The system is simulated for a period of 2 s. The shell is meshed with an 8 x 4 grid of 
quadratic shell elements, whereas the link is meshed with four cubic beam elements. The simulation was 
conducted using the ED scheme with p^ = 0, i.e. asymptotic annihilation, and a constant time step size 
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Fig. 6.7. Configuration of the snap-through problem. 



Fig. 6.8. System configurations at various instants in time. 






Fig. 6.9. Time history of the vertical displacements at points B (o), M (A) and E (+). ED scheme: solid line; EP 
scheme: dashed line . 
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Fig. 6.10. Time history of the bending moments Mn (°) and M 22 (A) at point E. ED scheme: solid line; EP scheme : 
dashed line . For clarity, the EP predictions were shifted down by 6 kN. 

A* = 2.5 1CT 3 S. 

During the first 90 degree rotation of the crank, the link pulls the shell downwards until snap-through 
takes place and curvature reverses. Curvature reversal initiates in the neighborhood of the link connection, 
then quickly propagates to the shell free edge, which undergoes violent oscillations. During the next 90 
degree rotation of the crank, the link now pushes the inverted shell upwards, until snap-through occurs 
and curvature reverts to its original sign. During the entire sequence, violent oscillations are observed, as 
depicted in fig. 6.8 that shows the system at various instants in time. 

Vertical displacements at point B, M, and E are shown in fig. 6.9. The displacement of point E closely 
follows the prescribed input imparted by the crank, whereas those at points B and M reflect the additional 
elastic vibrations of the shell. At first, the shell takes a double-S configuration: the curvature of the central 
part of the cylinder is already negative, while remaining positive along the simply supported edges. Snap- 



15 













































































Fig. 6.11. Time history of the lever mid-span axial force. ED scheme: top figure; EP scheme: bottom figure. Note the 
different scales for the two figures. 



Fig. 6.12. Time history of the lever mid-span bending moment. ED scheme: top figure; EP scheme: bottom figure. Note 
the different scales for the two figures. 


through occurs at about t & 0.7 s and vibrations occur in the inverted configuration. As the crank pushes 
the shell back up, local deformations appear at first, followed by a rapid snap-back to the original curvature 
at t « 1.25 s, when the crank is about to stop. Due to the speed of snap-back, violent oscillations about 
the original shell configuration are observed in the latter part of the simulation. The components Mu and 
M 22 of bending moment at point E are shown in fig. 6.10. Here again, the snap-through events are clearly 
identifiable. The lever mid-span axial force and bending moment are sown in figs. 6.11 and 6.12, respectively, 
Next, the system was simulated with the EP scheme, i.e . with poo = 1. In order to achieve convergence of 
the iterative solution procedure for the nonlinear equations of motion, the time step size had to be reduced 
to At = 1.25 HT 03 then 6.25 10“ 04 s at times t - 0.835 and 0.8475 s, respectively. A good agreement 
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configuration configuration 

Fig. 6.13. Configuration of the space antenna deployment problem. 

between the ED and EP predictions is shown in fig. 6.9, although discrepancies are apparent when significant 
oscillations occur. The bending moments shown in fig. 6.10 are once again in good agreement during the 
beginning of the simulation, but at the end of the snap-through phase ( i.e . t & 0.835 s), the EP scheme 
is unable to deal with the complex dynamic behavior of the system: a smaller time step is required for 
convergence, and violent high frequency oscillations of a purely numerical origin appear. The amplitudes 
of these oscillations are an order of magnitude larger than those predicted by the ED scheme. The same 
remarks can be made about the lever mid-span axial force and bending moments shown in figs. 6.11 and 6.12, 
respectively. In fact, the bending moment response, a simple superposition of oscillations involving the lowest 
two natural modes of the lever as predicted by the ED scheme, is completely obscured by numerical noise 
in the EP scheme. Since vibratory stresses are of great importance to designers, it is essential to assess 
the ability of new integration schemes to reliably predict these quantities. It is unfortunate that many 
scientific publications about geometric integration only present responses for preserved quantities such as 
total mechanical energy or momentum. The above plots demonstrate that while EP schemes might perform 
very well for the prediction of total energy, momentum, or even displacement fields, they are unable to 
reliably predict other important fields such as velocities and internal stresses. Consequently, such schemes 
are of little value in real life applications. 

6.3. Deployment of a Space Antenna. Consider a space antenna consisting of four square composite 
panels of size 5 x 5 m. A typical panel in shown in fig. 6.13. Each panel is reinforced along opposite edges 
AiBi and CiDi by beams of length L — 2 m. These reinforcing beams are attached at points and Bi to 
connector beams, and points C\ and D{ to re volute joints. The connector beams are attached to the re volute 
joints of the preceding panel. At points A\ and B \, Panel 1 is connected to the ground by means of revolute 
joints. The stowed configuration of the entire system is also depicted in fig. 6.13. 

The panels are made of laminated composite material with a 12 ply lay-up [0,90,45, -45,45, —45°]*, 
where the 90° direction is parallel to that of the reinforcing beams. The material properties of the composite 
are: longitudinal modulus El = 138 GPa, transverse modulus Et = 8.96 GPa, shearing modulus Glt = 
7.1 GPa, Poisson’s ratio vlt — 0.3, and ply thickness t p = 0.125 mm. The reinforcing and connector 
beams have a square cross-section of side a r = 10 mm and a c — 20 mm, respectively. They are made of a 
homogeneous isotropic material with the following properties: E — 140 GPa, and density p = 7000 kg/m 3 . 
Each of the eight revolute joints weighs 8 kg and is spring loaded to deploy the antenna. Each spring applies 
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Fig. 6.14. Configuration of the space antenna at various instants in time . 

Table 6.1 

Physical properties of the springs and dampers at the revolute joints. 


Revolute 

Joint at 

Mo 

[Nm] 

do 

[rad] 

n 

[Nms/rad] 

Ai,B\ 

1.0 

7t/2 

6 

20.0 

C 1 ,D 1 

0.8 

7T 

12 ' 

9.0 

c 2 ,d 2 

0.7 

7T 

12 

7.0 

Cs,D 3 

0.5 

7T 

12 

5.0 


a moment M s = Af 0 [ 1 - (0/0 o ) n ] at the joint. Furthermore, a damper is present in each joint and applies a 
moment Md — [i6, where 6 is the relative rotation at the joint. The constants Mo, 0o ? n an d A* are listed in 
table 6.1 for each joint. 

Under the effect of the springs in the revolute joint, the antenna deploys. Configurations of the system 
at various instants in time are shown in fig. 6.14, and the relative rotations at four revolute joints during 
deployment are shown in fig. 6.15. To validate the simulation, a simplified model of the system using 
geometrically exact beam elements to represent the panels was also run. The predictions of both models are 
in fair agreement. 

It should be noted that the composite lay-up used for the panels presents a 10% bending-twisting 
elastic coupling term. Consequently, the bending of the panels during deployment generates twisting, and 
the motion of the entire system becomes three-dimensional. To illustrate this effect, the rotations of the 
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Fig. 6.15. Relative rotation at the revolute joints. Joint at point B\: top-left, D\: top-right, D<i'- bottom-left, and D&: 
bottom-right figure. Shell model: solid line; beam model: dashed line . 






Fig. 6.16. Rotation of the reinforcing beams about the ii axis. Panel 1: top-left, Panel 2: top-right, Panel 3: bottom-left, 
Panel 4 •* bottom-right figure . 


reinforcing beams about the axis ii are shown in fig. 6.16. These rocking motions of up to 7 degrees generate 
large, mid-span shear forces in the connecting beams between Panels 1 and 2 that are shown in fig. 6.17. 
These effects are ignored in the beam model that predicts a two-dimensional motion. 

The bending moments M\\ at the center point of each panel are shown in fig. 6.18. Pronounced dis¬ 
crepancies are observed between the shell and beam models because the latter does not take into account 
the bending-twisting coupling behavior in the panels. Note the much higher frequency content predicted 
by the shell model. This is due to the twisting and transverse bending of the panels that occur at much 
higher frequencies than those associated with longitudinal bending. Significant twisting Myi and transverse 
bending M 22 moments develop in the panels, as depicted in fig. 6.19; of course, these twisting moments 
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Fig. 6.17. Transverse shear forces in the connector beams between Panels 1 and 2. Beam C 1 A 2 ’ solid line; beam D 1 B 2 : 
dashed line. 





Fig. 6.18. Mid-span bending moments Mu in the panels. Panel 1: top-left, Panel 2: top-right , Panel 3: bottom-left, 
Panel 4 ** bottom-right figure. Shell model: solid line ; beam model: dashed line . 


vanish in the beam model. Note that the transverse bending moment M 22 is of the same order of magnitude 
as the longitudinal bending moment Mu , confirming the plate-like nature of the deformation in each panel, 
an effect ignored by the beam model. 

7. Conclusions. In this work, a new geometric integration procedure was developed for the simulation 
of multibody system dynamics involving shells. The proposed scheme is “aware” of the qualitative features 
of the underlying partial differential equations. In particular, it evolves on the special manifold defined 
by two-parameter rotation fields, on the manifold defined by the presence of holonomic and non-holonomic 
constraints imposed by the mechanical joints and on the manifold of constant total linear and angular 
momenta. 

In contrast with the classical energy preserving approaches, the method presented here lets the system 


20 



























































































.!^. 

! 


: I : 




100 120 140 160 


TIME [sec] 




TIME [sec] 


Fig. 6.19. Mid-span twisting (Mi 2 ) and bending (M 22 ) moments in the panels. Panel 1: top-left, Panel 2: top-right , 
Panel 3: bottom-left, Panel 4 : bottom-right figure. 

drift away from the level of constant energy in a controlled manner. This achieves two important goals. 
First, a bound is placed on the total mechanical energy of the system at the discrete solution level, therefore 
achieving nonlinear unconditional stability. Second, the energy decay provides high frequency numerical 
dissipation. The proposed integrator can be used with arbitrary spatial discretizations, for example based 
on finite element or finite volume techniques. The shear locking phenomenon was controlled using the mixed 
interpolation of tensorial components approach. 

The proposed shell model was developed within the framework of a multibody dynamics analysis proce¬ 
dure that includes rigid bodies, beams, and various types of joints. The constraint forces are discretized so 
that the work they perform vanishes exactly at the discrete level. A shell revolute joint that connect shells 
to other elements of the model was developed. The efficiency and robustness of the proposed approach were 
demonstrated with specific numerical examples. 
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Appendix A. Rodrigues Parameters. 

A common representation of finite rotations [20] is in terms of Rodrigues parameters r = 2k tan<^/2, 
where (p is the magnitude of the finite rotation and k the components of the unit vector about which it takes 
place. The following notation is introduced tq = cos 2 <p/2 = 1 / (1 + r T r/4), and the finite rotation tensor 
R then writes 


^ Tq _ 

R(r) — I + ro r + — rr. 

The following decomposition of the rotation tensor is extensively used in this work 


R + I 


(A.l) 


(A.2) 


Appendix B. Orientation of a Unit Director. 

Consider a unit vector £ 3 , called a director , that rotates to a final orientation e 3 . For convenience, this 
director is considered to be the third unit vector of a triad S defined by £ x , £ 2 , hi rotating to a triad S* with 
orientation e 1? e 2 , e 3 . The relationship between these two triads is e a = R £ a , where R is an orthogonal 
rotation tensor. If one solely focuses on the director, this rotation tensor is not uniquely defined, as any 
rotation about the director leaves its orientation unchanged. A virtual change in the director orientation is 


Se 3 = e% Sip, 


where Sip is the virtual rotation vector, Sip — SRR T . 

The components of the virtual change in director orientation measured in S* become 


R T Se 3 = R T e%5'ip — 7^ R T 6ip — 1 % = 


sr 2 

sri , 

0 


(B.l) 


(B.2) 


where Sip* are the components of the virtual rotation vector in <S*. This relationship clearly demonstrates 
that arbitrary values of Sip 3 , corresponding to virtual rotations of the director about its own orientation, 
will not affect virtual changes in the director orientation, and hence, setting Sip 3 = 0 is a valid choice. The 
following notation is adopted 

Sip* = i^Sal +£2^2 = & M*; & = [ii,i 2 ]- (B. 3 ) 

Sa * is a 2 x 1 , “two parameter” virtual rotation vector. It follows that Sip_ — R Sip * = Rb Sa* , and hence 

Se a =Rl r a bSo?. (B.4) 


If Rodrigues parameters are used to parameterize R , an equivalent expression can be obtained for finite 
changes in director orientation with the help of eq. (A. 2 ) 

e af -e ai = R m 'iZbs* =Q m s*; L*=bs*, (B.5) 

where r* are the Rodrigues parameter measured in <S*, and s* the corresponding “two parameter” incremental 
rotation vector. 
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